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Abstract 

A procedure to obtain differentiation matrices is extended straightforwardly to yield new 
differentiation matrices useful to obtain derivatives of complex rational functions. Such 
matrices can be used to obtain numerical solutions of some singular differential problems 
defined in the complex domain. The potential use of these matrices is illustrated with the 
case of elliptic functions. 



1 Introduction 

In a sequel of papers and references therein), a Galerkin-type method has been 

used to solve boundary value problems and to find limit-cycles of nonautonomous dynam- 
ical systems. The method is based on the discretization of the differential problem by 
using differentiation matrices that are projections of the derivative in spaces of algebraic 
polynomials or trigonometric polynomials. This kind of matrices arise naturally in the 
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context of interpolation of functions and yield exact values for the derivative of polyno- 
mial functions at certain points selected as nodes. The class of functions that defines the 
domain of the differential operator determines the kind of differentiation matrix to be 
used. Thus, to get a discrete form of an operator acting on real functions that drop off 
rapidly to zero at large distances, we can use the skew-symmetric differentiation matrix 



constructed with the zeros xj of the Hermite polynomial Hj^{x) as lattice points^. On 
the other hand, to get a discrete form of an operator acting on periodic functions [Zj, the 
differentiation matrix 



should be used. In this case the lattice points can be chosen as the equidistant points 



Discrete forms of multidimensional differential operators can be obtained by using direct 
products of suitable one- dimensional differentiation matrices (see jH] for instance). 
Other approaches to obtain differentiation matrices can be found in [TU] . 
The purpose of this paper is to use the complex Hermite interpolation formula to find 
new differentiation matrices on the complex domain that produce an exact formula for the 
derivatives of complex rational functions and that can be used to approximate the solution 
of a singular differential equation in the complex domain. As examples, the derivatives 
of meromorphic two-periodic functions such as Jacobi elliptic functions and Weierstrass 
elliptic functions will be used to show the potential of this technique. 



2 Interpolatory differentiation matrices 

In order to illustrate the procedure to generate differentiation matrices, an already known 
matrix which gives the derivative of algebraic polynomials of the complex variable z 
is obtained by using the Hermite interpolation formula. 



"'^Equation gives an asymptotic expression for D_ 








= -TT + 2T13/N, J = 1, 2, ■ ■ ■ , 



N. 
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2.1 Algebraic polynomials 



Let f{z) be an analytic function on a domain G containing a closed rectifiable Jordan 
curve 7 and let Zk be N different points of 1(7) defining the polynomial 



N 



(^(z) = Y\iz - Zk). 



k=l 



Then, the unique polynomial p{z) interpolating f{z) associated to the set of points Zk is 
given by (0, i) 



(2.1) 



fiOuiO-uiz) 



2ni J Lu{() C - z 



The residual function R{z) = f{z) —p{z), z G 1(7), vanishes at Zk yielding that f{zk) 
p{zk)- To see this, write R{z) as 



(2.2) 



R{z) 



2m 



/(C) 



(C - ^)^(C) 



Since the integral of the right-hand side of this equation represents an analytic function 
on 1(7) we have that R{zk) = 0. 

To obtain the form of p{z), the integral in Eq. ()2.1|) can be calculated by the residue 
theorem. Since — u{z) is divisible hj ( — z, the integrand has simple poles only 
at zi,...,zn. Thus, the residue theorem yields the well-known Lagrange interpolation 
formula 



TV 



(2.3) 



n 



N 



Z - Zk 



Hk^ji^j - ^k) 



zeG. 



The differentiation matrix D associated to this formula can be obtained if the derivative 
of (12.311 is evaluated at Zj and written as 



(2.4) 
Since 



dp{zi 
dz 



N 



d 

dz 



' N N 
\\{Z-Zk)/ \\{Zj-Zk) 

■k^j k^j 



]\k^Mj-Zk) 

N 



' {Zi - Zk) ' 
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we get that the differentiation matrix for algebraic polynomials is given by 



( N 



1 



E 




z = 



J, 



(2.5) 



= < 




^ 7^ J- 



If f{z) is a polynomial of degree at most — 1, f{z) is given identically by p{z) and 
therefore formula ()2.4j) gives the exact derivative of f{z). Since f'{z) is another polynomial 
of this class, the derivatives of higher order can be obtained by applying successively the 
matrix D to the vector of values f{zj), i.e., 



is just the vector whose entries are f{zi). The functional form of f^'^\z) can be obtained 
through an interpolation of the values yielded by ()2.6|) . Since any set of different com- 
plex numbers belonging to G yields the same polynomial f{z), this result is independent 
of the points Zk- On the other hand, if f{z) is not a polynomial of degree at most N — 1 
a residual vector should be added to the right-hand side of ()2.4|) to get f'{zi). However, 
such a formula yields an good approximation to f'{zi) if the absolute value of the M-th 
term of the Taylor series of f{z) goes to zero rapidly as M is increased. 



2.2 Trigonometric polynomials 

The preceding arguments can be modified to consider the interpolation of periodic func- 
tions in terms of trigonometric polynomials. Let f{z) be a one-periodic analytic function 
with period 27i and let G be a domain of the open strip — vr < ^z < n, — oo < 'isz < oo, 
containing a closed rectifiable Jordan curve 7. 

Since any trigonometric polynomial t{z) = oq + J2T=i(^k cos kz + b/. sin kz) of degree at 
most m can be written in the form f{s) = s~''^q{s) under the the mapping s = ^p{z) = e*^ 
where q{s) is a polynomial of degree at most 2m in s, we need to take an odd number 
= 2m + l of different points Sk G (p{I{'j)), i.e, 2m + l different points Zk G /(7), to yield 
an exact interpolation formula in the case in which f{z) is a trigonometric polynomial of 
degree at most m. This fact can be shown as follows. 

Let us take = 2m + L The set of points Sk, k = 1,2, ... ,N define the polynomial 
^i^) = n^il-^ ~ Sk). The interpolant function p{s) to /(s) = f{ip~^{s)) corresponding 
to the set of A^ points Sk is given by 




(2.7) 
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where 7 = Since [s'^uj{()—('^uj{s)]/{( — s) is a polynomial in s of degree iV—1 = 2m, 

p{s) has the required form s~"^q{s), where q{s) is a polynomial of degree at most 2m, to 
represent a trigonometric polynomial t{z). 

To show that f{sk) = p{sk), let us consider the residual function R{s) = f{s) —p{s) which 
is now 



Ris) 



1 u{s) f /(C)C 



27rz J. (C - s)Cu{C) 



By definition, G = f{G) does not contains points Sfc(mod 27r) other than Sk therefore, 
the integral of the right-hand side of this equation represents an analytic function in 7(7) 
and we have that R{sk) = 0. 

Since s"^uj{() — ("^uj{s) is divisible hj ( — s, the poles of the integrand are simple and 
located at Sk- The residue theorem yields now 



N 



(2.8) 



q . \ m 



ULAs - Sk) 



j=l - Uk^jiSj - Sk) 



s e G 



and the trigonometric polynomial of degree m = (A^ — l)/2 that interpolates to f{z) is 



(2.9) 
Since 



p{z) 



N 



Zj}e 



i{N-l)(zi-z)/2 i-i-k^j 



n 



N 



J(N-l)(z,-z)/2 i-i-k^j 



n 



N 



p'^'^h 



nfc;,,(e-^-e-'=) 

. nf^,sin(^) 



zeG. 



n 



N 
k^j 



nM,sin(^) 



we obtain from (j2.9|) the Gauss interpolation formula 



N 



(2.10) 



p{z) 



A^ = 2m + l, zeG. 



By writing the derivative of this formula in the form given by ()2.4|1 we can obtain the 
differentiation matrix for trigonometric polynomials. Thus, the matrix D whose entries 
are given by 



(2.11) 



1 ^ 



cot 



1 

-CSC 



Zi Zj 



nN 



, Zj Zk • 



in terms of = 2m + 1 different points Zk G G, is a projection of d/dz in the subspace 
of trigonometric polynomials of degree at most (A^ — l)/2. Therefore, if f{z) belongs to 



5 



this space, p{z) = f{z), and f^'^\z) satisfies again an equation like ()2.6|) but in this case 
is the nth power of (EZH). The form of f^^'^z) can be obtained from the set of values 
f^"'\zi) through an interpolation. For a general one-periodic analytic function with period 
2tc a residual vector should be added to ()2.10|) . 

It is worth to notice that in the case in which the set of points Zk are real numbers, the 
matrix ()2.1H) becomes the matrix used previously in [H], [7j. However, if the points lay 
on a straight line which is parallel to the imaginary axis, D takes the form 



D.J 



N 

coth 



Vi - Vk 



-csch 

2 V 2 



nfe^jSmh(^— ) 
nfc^.smh(^- ) 



^ 7^ J, 



where yk = '^Zk, at the same time that the polynomial to differentiate becomes a linear 
combination of real hyperbolic functions. 



2.3 Rational functions 

Equation ()2.7|1 suggests the form of a interpolant rational function in the case in which 
f{z) is a meromorphic function. 

Let G be a domain that contains a closed rectifiable Jordan curve 7 and let f{z) be a 
meromorphic function with only one pole at 2; = a ^ G of order m. Let us choose N 
different points z^ of 1(7) and construct the polynomial u^z) = Y[k=ii^ ~ ^k)- Thus, the 
rational function interpolating to f{z) corresponding to the set of points Zk is given by 

2.12 P{z) = — / — — dC. 

This can be shown by the same arguments used in the previous case. Again, {z—a)^uj{Q — 
(C — a)'^uj{z) is divisible hy C, — z. Let K^{z, () be such a quotient, i.e., 

Kn[z,Q = • 

Since KYf{z,() is a polynomial of degree — 1 in z, p{z) is a rational function of form 
q{z)/{z — a)"^ that interpolates to f{z) at Zk, where q{z) is a polynomial of degree at most 
— 1. The residual function 

R(z) = [ IM^^dC 
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vanish at z^. because the integral is an analytic function in 1(7) and a ^ G. Therefore, 
p{zk) = fi^k)- The residue theorem yields now 



N 



(2.13) 



z — a 



Uk^jiZj - Zk)' 



zeG. 



The derivative of this equation at Zi can be written in the form ()2.4j) where we have now 



' ^ 1 



(2.14) 



m 



k^i 



V 



iz,-ar/iz,-ar u;'{z,) 
Zi-Zj uj'{zj)'- 



« =J, 



i 7^ J- 



Obviously, if f{z) is a rational function of the form q{z)/{z — a)™ where q{z) is a poly- 
nomial of degree at most — 1, f{z) = p{z) and formula (j2.4p becomes 



(2.15) 



dfiz,i) 
dz 



N 



However, the powers of D do not give the derivatives of higher order as in the previous 
cases since f'{z) does not has the required form: it has a pole of order m + 1 at z = a. 
Despite this, it is possible to obtain the value of f^"\z) at Zk by using the matrix 



(2.16) 



DJm) = Dr 



m+n—lDfii+n—2 



Each matrix is defined by (|2.14p where the parameter m, defining the order of the 
pole, is substituted by each value of the index k. Thus, it should be clear that if f{z) has 
the form given above, ()2.6|) becomes 

(2.17) /(") = D„(m)/, n = 0,1,2,.... 

This equation can be generalized to the case in which f{z) is a rational function of form 

Z^fc=0 ^kZ 



(2.18) 



f{z) 



[z - ai)™i {z - a2)"2 ■■■{z- arY 



where ^ G, / = 1, 2, ■ ■ ■ , r. A straightforward calculation gives the generalized form of 

(2.19) =Z}„(mi,m2,...,m,)/, n = 0,1,2,..., 

where now Dn{mi, m2, . . . ,mr) stands for the ordered matrix product 



(2.20) 



Dn{mi,m2,. . . ,mr) = JJ Dmi+k,m2+k,...,mr+k 



k=n—l 



and D 



is the matrix whose entries are given by 



5Z _ 



(2.21) 



n 



Zi-Zj uj'[Zj) \Zi-ak 



^ = J, 



It should be noticed that ()2.19p is an exact formula whenever N > M + nr — 1, where 
M is the degree of the polynomial in the numerator of ()2.18p . The reason is that after 
each differentiation the numerator of the derivatives of f{z) is a polynomial whose degree 
grows by r. If the function f{z) to differentiate has poles at ai, . . . , of orders ni, . . . , 



, rrir 



with Hi < mi, 



Tir < rrir 



formula ()2.19|) is still exact whenever 



instead mi, 

> M + nr — 1 + X]I=i("^fc ~ ^k)- The reason for this is that in such a case f{z) 
can be writen in the form ()2.18|) where the numerator is now a polynomial of degree 
^ + Sfc=i("^fc ~ "^fc)- Obviously, in this case is much better to use the differentiation 



matrix D„(ni, n2, 



instead D„(mi,m2, 



, rrir 



for numerical purposes. 



3 Applications 

As stated before, the numerical solution of differential problems can be accomplished by 
the use of differentiation matrices, and in the case of a differential problem in the complex 
domain, the differentiation matrices introduced in this paper may be useful. To illustrate 
the potential of their use, we choose two meromorphic cases which are important in ap- 
plications: Jacobi elliptic functions and Weierstrass elliptic functions. In both cases it is 
possible to establish the numerical convergence of the results since the derivatives of these 
functions are known. We also obtain approximate solutions of a differential equation with 
a regular singularity at z = 0. Before beginning these examples is convenient to alert the 
reader to the fact that the numerical implementation of the matrices for rational functions 
given above may need a high-precision code: in most cases, the usual 16-digit precision is 
not enough to obtain accurate results. 



3.1 A rational function 

Let us consider the function 

(3.1) f{z) 



Z^ + Z + 1 
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Table 1: The norm of the relative error = \{f"' — p'")/ f"'\oo for 
in terms of the number of nodes A^. 





En 


4 


0.657 


5 


0T36 


6 


0.155 xlO-i 


7 


0.742 xlO-3 


8 


0. xlQ-^^ 


9 


0. xlO-i*^ 


10 


0. xlO~^^ 


11 


0. xlQ-^^ 



According to the results of the last section, to obtain the exact value of the n-th derivative 
of (j3.H) . [cf. Eq. (j2.17p ]. it is necessary to choose N > 7 different points 7^ in the 
complex plane to build the matrix ()2.16|) where obviously, a = and m = 10. 
As an example, let us take the third derivative of ()3.1|) . Table ITT] shows the max-norm of 
the vector whose entries are \ [f^^\zj) — '^^^-^^{D3{10))jkf{zk)]/f^^\zj)\ (the relative error) 
in terms of A^. The nodes were chosen to be as evenly spaced on the ray z = {1 + i)t, 
1/2 < t < 1, i.e., Zk = + k/N)/2, k = 1, . . . , N, and the computations were made 

with the standard 16-digit precision. As it can be seen from Table ITTl the error vanishes 
for values of A^ greater than 7, as expected. 



3.2 Elliptic functions 

As is known, an elliptic function is a doubly periodic function which is analytic except at 
poles [21-CO of the two simplest cases of elliptic functions corresponds to Jacobi's 

functions; the other corresponds to Weierstrass' P-function. 

The differentiation matrices given above do not apply to functions like these, however it 
is possible to build a matrix which is expected to provide approximate values of f^^\zk) 
along straight lines inside the fundamental paralelograms and near to the poles. Since an 
elliptic function f{z) becomes a one-periodic function if z is constrained to move along 
a straight line defined by one of the periods, such a matrix can be constructed by using 
trigonometric polynomials divided by an algebraic polynomial with zeros of suitable orders 
taken at the poles of f{z). This procedure is equivalent to taking apart the matrix ()2.21|) 
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and incorporating only the singular terms in ()2.11|) to yield the new matrices 



(3.2) Dn{mi,m2,...,mr)= JJ 

k=n—l 

where 



(3.3) (A 




To test the numerical performance of this matrix we choose two numerical examples. The 
first one corresponds to Jacobi's function f{z) = sn(z| | ) which has two periods AK and 
2iK' and two simple poles at iK' and 2K + iK', where 

K = / — . K' = K^ 1.854. 

io Vl-(sin'^)/2 

Therefore, to build the differentiation matrix we take 27r-periodic trigonometric polyno- 
mials, r = 2, /ii = /i2 = 1, tti = iK' , and a2 = 2K + iK' in (j3.3p . and to measure the 
approximation of f'{zj) by p'{zj) = J2k=ii^^i^^ ^))jkf{zk) we use the max norm. Here, 
f'{z) = cia{z\ I ) dn(2;| |). The results are displayed in Figure 1, where the max-norm 
of the error f'{zj) — p'{zj) is plotted against the number of nodes showing numerical 
convergence. The number of digits of precision used in the calculations is 16 and the 
nodes are chosen to be as evenly spaced on the ray z = {2 + i)t, 1/2 < t < 1, i.e., 
Zk = (2 + i){l + k/N)/2, k = l,...,N. The norm of the error is 1.6 x 10^^ for = 10 
and 1.0 X 10"^ for N = 20. 

Our second example is the Weierstrass P-function, which has a double pole at the origin 
and two periods coi and uj2. The derivative of V{z) satisfies 

(P'(^))2 = 4(P(^))3_^^p(^)_^3, 

where g2 and gs are the elliptic invariants which are related to both uJi and uj2 jlj- 
To approximate the first derivative of V{z) at the nodes, we need to take n = 1 and r = 1 
in ()3.2j] and fii = 2 and ai = in ()3.3j) . The precision of the calculations, the nodes and 
the kind of trigonometric polynomials used to construct the differentiation matrix are the 
same as above. The numerical results are displayed in Figure 1 and show again numerical 
convergence with a small number of nodes. For = 10 the norm of the error is 1.5 x 10~^ 
and 1.0 X 10"^ for = 20. 
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Figure 1: The norm of the error Ei\i = \ f' — p'\oo against the number of nodes for 
the Jacobian ehiptic function (soUd hne) and Weierstrass' P-function (broken hne). 
The vertical axes are scaled accordingly to Jacobi's case (the left-hand side axis) or 
Weierstrass' case (the right-hand side axis). 



3.3 Krummer's equation 

As a final example, we consider the Krummer differential equation, written as an eigen- 
value problem 

(3.4) ^^ + (^-^)^ = «^^^)' 

which has a regular singularity at z = and an irregular singularity at oo. As is well 
known, the single- valued solution of this equation is the confluent Hypergeometric function 
M(a, 6, z) =iFi(a, 6, z). Since this function can be approximated by algebraic polynomials 
for h 7^ —n in a positive integer), we can obtain approximate solutions of this differen- 
tial equation by using the differentiation matrix D given by (|2.1H) to approximate the 
derivative of M(a, 6, z) and solving the A^-dimensional eigenvalue problem 

(3.5) L/a = A/a, L = + (^^ - Z)D 

where Z is a. diagonal matrix whose nonzero elements are the nodes zi, . . . , zn, & is a 
complex number {b ^ —n), Iat is the identity matrix of dimension and the eigenvalue 
A is the value of a at which M(a, b, Zk) is to be approximated by (/A)fc- Let us denote by 
Ma the vector whose kth component is M{a, b, Zk). Since the nth coefficient of the power 
series of M(a, 6, z) is given by 

a(a + l)(a + 2) . . . (a + n - 1) 



b{b + l){b + 2) ...{b + n-l)n\' 
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the best approximation obtained for a given set of parameters b, N, Zk, is given by the 
eigenvector fx (A = a) corresponding to the eigenvalue with lowest absolute value A^- To 
construct the matrix L in (j3.5p we choose = 21 and = 5(1 + i)k/N, k = 1, . . . N. 
For b we take the cases 6 = 5/2 and b = 3 + 2i. In order to compare the approximate 
and exact results, we normalize both vectors Mxm and fxm with the max-norm. The 
calculations were made with 16 digits of precision and the results are displayed in Fig. 2. 
The absolute errror \\Mxm - /Am||o is 0.0675659 for b = 5/2 and 0.0426948 for 6 = 3 + 2i. 




Figure 2: Normalized real (Re) and imaginary (Im) parts of fxm plotted versus their 
index (solid lines). Case (a) corresponds to 6 = 5/2 and = —0.301513 + 1.00758i 
and case (b) to 6 = 3 + 2i and Am = —0.381925 + 0.5275331 They are compared 
with the exact values of the Krummer function (broken lines). The matrix L is 
constructed with 21 nodes Zk = 5(1 + i)k/N, k = 1, . . . 21. The scale on the left- 
hand side vertical axis corresponds to the real part and the one on the right-hand 
side to the imaginary part. 



4 Concluding remark 



According to the results of section |2l the process of interpolation in vector spaces of 
polynomials of dimension maps the derivative d/dx into a. N x N matrix D. The fact 
that a differential operator acting on a vector space of finite dimension can be written 
as a matrix is not a surprise of course, however it should be noted that the matrix D 
yields the derivative of a function by taking the values of the function at A^ arbitrary (but 
different) points including the point where the derivative is to be evaluated, i.e., it acts on 
a function as a nonlocal operator in spite of the local character of a differential operator 
as the derivative. 
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